# This script analyzes the results of LC_processVS_hpc.R

Packages <- c("rstan", "coda", "ggmcmc", "bayesplot", "tidyverse", 
              "loo", "purrr")
suppressMessages(invisible(lapply(Packages, library, character.only=TRUE)))
theme_set(theme_bw())

setwd("C:/Users/tms1044/Desktop/lc_mod")

n.cores <- 4
chn.dir <- "out/"
sum.dir <- "summaries/"
mods <- c("Cens", "Clim", "ClimCens", "Topo", 
          "TopoCens", "TopoClim", "TopoClimCens")

# GGMCMC
ac_cols <- c(rgb(0,0,0,0.5), rgb(1,0,0,0.5), rgb(0,1,0,0.5), 
             rgb(0,0,1,0.5), rgb(1,1,0,0.5), rgb(1,0,1,0.5))
beta.fls <- list.files(sum.dir, pattern="out_betas_", full.names=TRUE)
names(beta.fls) <- mods
map(beta.fls, readRDS) %>%
  map(., ~(ggs_density(ggs(.x)) + facet_wrap(~Parameter, scales="free")))
## $Cens

## 
## $Clim

## 
## $ClimCens

## 
## $Topo

## 
## $TopoCens

## 
## $TopoClim

## 
## $TopoClimCens

map(beta.fls, readRDS) %>%
  map(., ~(ggs_traceplot(ggs(.x)) + facet_wrap(~Parameter, scales="free")))
## $Cens

## 
## $Clim

## 
## $ClimCens

## 
## $Topo

## 
## $TopoCens

## 
## $TopoClim

## 
## $TopoClimCens

map(beta.fls, readRDS) %>% map(., ~(ggs_autocorrelation(ggs(.x)) + 
                                      facet_wrap(~Parameter) +
                                      scale_fill_manual(values=rep(NA,6)) +
                                      scale_colour_manual(values=ac_cols)))
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## $Cens

## 
## $Clim

## 
## $ClimCens

## 
## $Topo

## 
## $TopoCens

## 
## $TopoClim

## 
## $TopoClimCens

map(beta.fls, readRDS) %>% map(., ~(ggs_geweke(ggs(.x))))
## $Cens

## 
## $Clim

## 
## $ClimCens

## 
## $Topo

## 
## $TopoCens

## 
## $TopoClim

## 
## $TopoClimCens

map(beta.fls, readRDS) %>% map(., ~(ggs_Rhat(ggs(.x))))
## $Cens

## 
## $Clim

## 
## $ClimCens

## 
## $Topo

## 
## $TopoCens

## 
## $TopoClim

## 
## $TopoClimCens

map(beta.fls, readRDS) %>% 
  map(., ~(ggs_pairs(ggs(.x), lower=list(continuous="density"))))
## $Cens

## 
## $Clim

## 
## $ClimCens

## 
## $Topo

## 
## $TopoCens

## 
## $TopoClim

## 
## $TopoClimCens

map(beta.fls, readRDS) %>%
  map(., ~(ggs_running(ggs(.x)) + facet_wrap(~Parameter, scales="free")))
## $Cens

## 
## $Clim

## 
## $ClimCens

## 
## $Topo

## 
## $TopoCens

## 
## $TopoClim

## 
## $TopoClimCens

map(beta.fls, readRDS) %>%
  map(., ~(ggs_compare_partial(ggs(.x)) + facet_wrap(~Parameter, scales="free")))
## $Cens

## 
## $Clim

## 
## $ClimCens

## 
## $Topo

## 
## $TopoCens

## 
## $TopoClim

## 
## $TopoClimCens

# effective sample size
par(mfrow=c(3,3))
map(beta.fls, readRDS) %>% 
  walk(., ~hist(effectiveSize(.), main="", xlab="Effective Sample Size"))
par(mfrow=c(3,3))

map(beta.fls, readRDS) %>% 
  walk(., ~hist(effectiveSize(.)/(250*6), main="",
                xlab="n_eff/n_iter", xlim=c(0,1.5)))
par(mfrow=c(1,1))

# WAIC
waic.ls <- map(list.files(sum.dir, pattern="waic", full.names=TRUE), readRDS) 
names(waic.ls) <- mods
map_dbl(waic.ls, ~(.$waic)) %>% sort
##         Topo     TopoCens TopoClimCens     TopoClim         Clim 
##    -796749.0    -795971.0    -785065.4    -772394.8    -767220.5 
##     ClimCens         Cens 
##    -763440.7    -759462.5
compare(waic.ls[[1]], waic.ls[[2]], waic.ls[[3]], waic.ls[[4]], 
        waic.ls[[5]], waic.ls[[6]], waic.ls[[7]])
##              waic      se_waic   elpd_waic se_elpd_waic p_waic   
## waic.ls[[4]] -796749.0    1075.3  398374.5     537.6      84739.9
## waic.ls[[5]] -795971.0    1105.9  397985.5     552.9      84307.0
## waic.ls[[7]] -785065.4    1116.2  392532.7     558.1      77443.9
## waic.ls[[6]] -772394.8    1087.3  386197.4     543.7      80668.8
## waic.ls[[2]] -767220.5    1082.6  383610.2     541.3      83730.8
## waic.ls[[3]] -763440.7    1110.2  381720.4     555.1      84788.8
## waic.ls[[1]] -759462.5    1104.8  379731.2     552.4      81700.1
##              se_p_waic
## waic.ls[[4]]      24.9
## waic.ls[[5]]      28.6
## waic.ls[[7]]      26.5
## waic.ls[[6]]      24.3
## waic.ls[[2]]      24.9
## waic.ls[[3]]      27.3
## waic.ls[[1]]      29.6
# LOO
loo.ls <- map(list.files(sum.dir, pattern="loo", full.names=TRUE), readRDS) 
names(loo.ls) <- mods
map_dbl(loo.ls, ~(.$looic)) %>% sort
##         Topo     TopoCens TopoClimCens     TopoClim         Clim 
##    -753354.2    -751749.9    -740253.9    -728311.2    -723893.4 
##     ClimCens         Cens 
##    -720686.0    -715829.4
compare(loo.ls[[1]], loo.ls[[2]], loo.ls[[3]], loo.ls[[4]], 
        loo.ls[[5]], loo.ls[[6]], loo.ls[[7]])
##             looic     se_looic  elpd_loo  se_elpd_loo p_loo     se_p_loo 
## loo.ls[[4]] -753354.2    1080.4  376677.1     540.2    106437.3      57.2
## loo.ls[[5]] -751749.9    1110.6  375875.0     555.3    106417.5      58.0
## loo.ls[[7]] -740253.9    1120.6  370126.9     560.3     99849.7      57.4
## loo.ls[[6]] -728311.2    1091.6  364155.6     545.8    102710.6      56.6
## loo.ls[[2]] -723893.4    1086.5  361946.7     543.3    105394.3      57.5
## loo.ls[[3]] -720686.0    1115.1  360343.0     557.5    106166.2      58.1
## loo.ls[[1]] -715829.4    1108.8  357914.7     554.4    103516.6      58.1
map(loo.ls, plot)

## $Cens
## NULL
## 
## $Clim
## NULL
## 
## $ClimCens
## NULL
## 
## $Topo
## NULL
## 
## $TopoCens
## NULL
## 
## $TopoClim
## NULL
## 
## $TopoClimCens
## NULL
# OOS prediction
mspe.ls <- map(list.files(sum.dir, pattern="MSPE", full.names=TRUE), readRDS)
names(mspe.ls) <- mods
map_dbl(mspe.ls, sqrt) %>% sort
##         Clim     ClimCens         Cens     TopoClim TopoClimCens 
##    0.1810194    0.1811403    0.1854589    0.1893665    0.1899082 
##         Topo     TopoCens 
##    0.1909966    0.1914862
lpd.ls <- map(list.files(sum.dir, pattern="lpd", full.names=TRUE), readRDS)
names(lpd.ls) <- mods
map_dbl(lpd.ls, c) %>% sort
##         Clim     ClimCens         Topo         Cens     TopoCens 
##    -2.170384    -2.139518    -2.103142    -2.101887    -2.082612 
##     TopoClim TopoClimCens 
##    -2.031828    -2.006136
# Gelman
par(mfrow=c(3,3))
map(list.files(sum.dir, pattern="gelman", full.names=TRUE), readRDS) %>% 
  walk2(.x=., .y=mods, ~({hist(x=.$psrf[,2], main=.y, xlim=c(0.95, 1.2)); 
    abline(v=1.1, lty=3, col="red")}))
map(list.files(sum.dir, pattern="out_betas", full.names=TRUE), readRDS) %>% 
  walk2(.x=., .y=mods,  ~(gelman.plot(x=.x, ylab=paste(.y, "shrink factor"),
                                      ylim=c(0.95, 1.2))))

# Geweke
gew.ls <- map(list.files(sum.dir, pattern="geweke", full.names=TRUE), readRDS) 
for(m in 1:length(mods)) {
  par(mfrow=c(2,3))
  for(i in 1:length(gew.ls[[m]])) {
    plot(density(gew.ls[[m]][[i]]$z), main=mods[m], xlim=c(-5,5), ylim=c(0,0.5))
    curve(dnorm(x, 0, 1), from=-5, to=5, add=TRUE, col="dodgerblue")
    abline(v=c(-1.96, 1.96), lty=3)
  }
}

map(list.files(sum.dir, pattern="out_betas", full.names=TRUE),
    readRDS) %>% walk(geweke.plot)